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Abstract —  Signal  detection  in  noisy  data  set  is  a  common 
problem  in  signal  processing.  Detection  of  the  hemody¬ 
namic  response  function  (HRF)  embedded  in  randomized 
event-related  fMRI  (rER-fMRI)  time  series  is  an  example 
of  this  problem.  So  far,  most  studies  that  set  out  to  obtain 
unbiased  HRF  use  some  forms  of  time-window  (TW)  av¬ 
eraging  method  to  extract  HRF  from  the  rER-fMRI  data. 
In  this  paper  we  applied  two  methods,  Cepstral  Analysis 
and  Conjugate  Gradients  (CG)  for  Deconvolution  to  extract 
HRF.  These  methods  depend  only  on  the  knowledge  of  when 
events  occured  and  do  not  require  any  apriori  information 
about  the  HRF.  These  methods  and  the  popular  TW  av¬ 
eraging  method  are  tested  on  simulated  data,  as  well  as  in 
vivo  data  obtained  from  rER-fMRI  experiments.  All  three 
methods  identified  timing  of  HRF  accurately,  but  only  CG 
for  Deconvolution  method  was  robust  in  reproducing  the 
shape  under  varying  experimental  conditions. 

Keywords — fMRI;  Event-related;  Hemodynamic  response; 
MRI;  Conjugate  Gradients;  Cepstrum;  ISI. 


I.  Introduction 

The  randomized  event-related  fMRI  (rER.-fMRI)  experi¬ 
mental  designs  have  become  increasingly  popular  in  recent 
years  as  the  designs  have  drastically  improved  efficiency 
[1]  over  the  fixed  inter-stimulus  interval  (ISI)  designs.  The 
rER.-fMRI  designs  allow  different  stimuli  to  be  presented 
in  random  sequences  and  timing,  and  thus,  eliminating  po¬ 
tential  confounds  caused  by  anticipation,  habituation,  or 
other  strategy  effects  [2]. 

The  fMRI  signal  is  caused  by  changes  in  blood  oxygena¬ 
tion  and  blood  volume  that  result  from  the  neural  activity 
induced  by  the  stimuli.  The  timing  characteristics  of  this 
vascular  response,  or  the  hemodynamic  response  function 
(HRF),  can  be  extracted  from  the  rER.-fMRI  data  assum¬ 
ing  linearity  between  the  fMRI  signal  and  the  underlying 
neuronal  activity.  The  linearity  has  been  demonstrated  by 
various  studies  [3]- [5]  in  the  intermediate  ISI  range  (2-15s), 
but  not  always  true  for  short  ISI  [6]. 

So  far,  most  studies  that  set  out  to  obtain  unbiased  HRF 
use  some  forms  of  time-window  (TW)  averaging  method  to 
extract  the  HRF  from  the  rER.-fMRI  data.  We  found  that 
the  TW  averaging  method  could  be  unreliable,  depending 
on  the  mean  ISI  (mISI)  and  randomness.  The  goal  of  this 
study  was  to  develop  robust  algorithms  to  extract  HRF 
with  better  accuracy  than  the  TW  averaging  method. 

II.  Materials  and  Methods 

A.  Theories 

By  assuming  a  linear  time-invariant  model  for  the  ob¬ 
served  fMRI  response,  rER.-fMRI  signal  y(t)  can  be  mod¬ 
eled  as  the  convolution  y{t)  =  s(t)  *  h{t)  of  the  HRF  h{t) 
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and  event  onset  time  series  s(<),  which  is  represented  as  a 
series  of  impulses,  s(t)  =  [10010..]. 

A.l  Time  Window  Averaging 

TW  averaging  is  a  common  method  used  to  extract  the 
HRF  [7].  HRF  is  estimated  by  using  a  fixed  time  window 
starting  with  each  stimulus  and  then  averaging  the  signals 
measured  within  each  of  these  time  windows.  The  idea  is 
to  increase  low  signal  to  noise  ratio  of  rER.-fMRI  signals  by 
averaging. 


A. 2  Cepstral  Analysis 

Cepstral  analysis,  a  nonlinear  signal  processing  tech¬ 
nique,  is  commonly  applied  to  signals  that  have  been  com¬ 
bined  by  convolution,  such  as  in  speech  processing  and  ho¬ 
momorphic  filtering  [8].  Complex  cepstrum  (CCeps)  of  x(t) 
can  be  calculated  by  taking  logarithm  of  the  FFT  of  x(t) 
and  applying  inverse  FFT  to  the  result.  The  transforma¬ 
tion  of  a  signal  into  its  CCeps  enables  us  to  use  an  operation 
that  has  the  same  algebraic  properties  as  addition,  instead 
of  using  multiplication  in  the  frequency  domain. 

We  were  able  to  extract  h(t)  via  cepstral  analysis  by 
subtracting  the  CCeps  of  s(t)  from  the  CCeps  of  y(t.)  and 
taking  the  inverse  CCeps  of  the  result.  Then  a  low-pass 
filtering  was  applied  to  obtain  a  smooth  HRF. 


y(t) 
Y(f) 
log(Y(f)) 
F~1{log(H(f))} 
CCeps[h{t )] 
h(t) 


s(t)  *  h(t)  TimeDomain 
S(f).H(f)  Frequency  Domain 

log(S(f))+log(H(f)) 
F-1{«ofl(F(/))}-F-1{%(5(/))} 
CCeps[y(t )]  —  CCeps[s(t)\ 
CCeps~1{CCeps[y(t)\  —  CCeps[s(t )]} 


A. 3  Conjugate  Gradients  (CG)  for  Deconvolution 

The  CG  method  is  an  iterative  method  for  solving  a  lin¬ 
ear  system  of  equations  Ax  =  6,  with  symmetric  positive 
definite  (SPD)  coefficient  matrix  A.  However,  with  simple 
modifications  this  method  can  be  applied  to  nonsymmetric 
and  overdetermined  least-squares  problems  [9]. 

Following  the  previous  notation,  convolution  can  be  ex¬ 
pressed  as  Y (/)  =  S(f).H(f)  in  the  Fourier  domain.  In  or¬ 
der  to  solve  the  deconvolution  problem  of  extracting  HRF, 
we  applied  the  CG  algorithm  [9],  [10]  to  compute  a  least- 
squares  solution,  (see  the  detailed  algorithm  in  the  Ap¬ 
pendix) 

min  ||  Y(f)  —  ||2 
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Fig.  1.  Simulated  rER-fMRI  signal  with  5s  mISI,  generated  by  con¬ 
volving  the  ideal  HRF  with  the  event  onset  time  series  and  sam¬ 
pling  at  2s  interval. 

B.  Simulations 

An  ideal  HRF  [11]  was  used  as  the  impulse  response  func¬ 
tion  of  the  model.  The  stimulus  onset  time  series  were  gen¬ 
erated  with  randomized  ISI  for  a  given  mISI.  Then  the  ideal 
rER.-fMRI  signal  was  simulated  by  convolving  the  ideal 
HRF  and  the  stimulus  onset  time  series  (Fig  .1).  Finally 
the  ideal  rER-fMRI  signal  was  sampled  at  an  interval  of 
2s,  the  same  as  the  repetition  time  (TR)  used  in  the  fMRI 
experiments.  All  simulation  data  were  kept  within  a  total 
time  series  of  240s. 

Two  types  of  simulations  were  performed  to  compare  the 
three  HRF  extraction  algorithms.  First,  simulated  data 
with  different  mISI  values:  2s,  5s  and  10s,  were  used  to 
test  the  sensitivity  of  the  HRF  extraction  algorithms  to  a 
range  of  mISI.  Secondly,  simulated  data  with  5s  mISI  but 
different  randomness  were  used  to  test  the  stability  of  the 
three  methods. 

C.  rER.-fMRI  experiments 

The  five  stimulus  onset  time  series  (2s,  5s,  and  20s  mISI, 
and  the  additional  two  5s  mISI  series  with  different  ran¬ 
domness)  implemented  for  the  simulation  data  were  used  to 
construct  visual  stimulation  paradigms  for  the  rER.-fMRI 
experiments.  Each  stimulus  consists  of  a  black-and-white 
checkerboard  displayed  for  250ms  at  every  event  onset  time. 
The  experiments  were  performed  on  a  clinical  1.5T  MRI 
scanner  by  using  a  T2*  gradient-echo  2D  spiral  sequence 
with  TR  of  2s.  Twenty  slices  were  acquired  per  TR  to 
cover  the  whole  brain.  The  spatial  resolution  was  3.75mm 
X  3.75mm  X  7mm. 

The  fMRI  data  were  processed  and  analyzed  by  using 
Statistical  Parametric  Mapping  [11].  The  image  volumes 
were  adjusted  for  slice-timing,  realignment,  normalization, 
and  smoothing  prior  to  statistical  analysis  for  the  signifi¬ 
cance  of  activation.  In  order  to  obtain  meaningful  fMRI 


mean-ISI  =  10  sec  Estimated  HRF 


Fig.  2.  (Simulation)  The  effect  of  decreasing  mISI  on  HRF  extraction 
algorithms.  Simulated  signal  time  courses  are  shown  on  the  left 
column.  Results  of  the  HRF  extraction  algorithms  are  plotted  on 
the  right  column.  The  ideal  HRF  (dot-dashed)  is  reproduced  very 
well  by  the  Cepstral  Analysis  (dashed)  and  the  CG  for  Deconvo¬ 
lution  (solid)  methods,  but  not  very  well  by  the  TW  Averaging 
(dotted)  method. 

response  for  the  comparison  of  the  algorithms,  three  adja¬ 
cent  voxels  with  significant  activation  (t-score  >  4)  were 
identified  at  the  same  location  in  the  primary  visual  cortex 
from  all  data  sets.  The  raw  time  course  was  temporally 
smoothed  by  a  zero-mean  Gaussian  with  2s  standard  devi¬ 
ation.  The  average  time  course  of  the  three  voxels  for  each 
data  set  was  used  to  extract  HRF  with  the  three  proposed 
methods.  Note  that  these  methods  were  also  robust  for  a 
single  voxel  time-course  analysis. 

III.  Results 

A.  Simulations 

Cepstral  Analysis  and  CG  for  Deconvolution  methods 
worked  almost  perfectly  for  the  zero-noise,  simulated  rER.- 
fMRI  signals  for  mISI  of  2s,  5s,  and  10s  (Fig.  2).  Maximum 
mean-square  errors  (MSE)  were  calculated  as  0.15  and  0.09 
for  the  Cepstral  Analysis  and  CG  for  Deconvolution  meth¬ 
ods,  respectively.  On  the  other  hand,  results  of  TW  Aver¬ 
aging  method  were  unstable.  MSE  was  calculated  as  high 
as  20.51  for  5s  mISI.  Also  undershooting  or  overshooting 
of  the  peak  was  observed  for  all  cases. 

Cepstral  Analysis  and  CG  for  Deconvolution  methods 
were  very  robust  for  5s  mISI  of  various  randomness  (Fig. 
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mean-ISI  =  5  sec  &  Type  -  B  Estimated  HRF  mean-ISI  =  5  sec 


Time  (sec)  Time  (sec)  Time  (sec) 


mean-ISI  =  5  sec  &  Type  -  C  Estimated  HRF  mean-ISI  =  2  sec 


Time  (sec)  Time  (sec)  Time  (sec) 


Fig.  3.  (Simulation)  The  effect  of  changing  randomness  on  HRF 
extraction  algorithms.  Simulated  signal  time  courses  for  5s  mISI 
are  shown  on  the  left  column.  Results  of  the  HRF  extraction 
algorithms  are  plotted  on  the  right  column  with  the  same  line 
style  notations  as  those  in  Fig. 2.  The  TW  Averaging  method 
produces  the  least  accurate  results  among  the  three  methods. 


Fig.  4.  (In  vivo  data)  The  effect  of  decreasing  mISI  on  HRF  extrac¬ 
tion  algorithms.  The  in  vivo  time  courses  are  shown  on  the  left 
column.  Results  of  the  HRF  extraction  algorithms  are  plotted  on 
the  right  column  with  the  same  notations  as  those  in  Fig. 2.  Note 
that  the  result  of  Cepstral  Analysis  for  10s  mISI  was  cut  at  -0.7, 
although  it  has  a  minimum  value  of  -1.1 


3)  with  maximum  MSE  values  0.15  and  0.09,  respectively. 
But  the  performance  of  the  TW  Averaging  method  was  un¬ 
stable  even  for  the  same  mISI.  From  simulations,  we  found 
that  TW  Averaging  performs  well  in  estimating  time-to- 
peak  (TTP)  of  HRF,  but  the  method  is  not  reliable  in  re¬ 
producing  the  shape,  such  as  maximum  percentage  change 
(MPC)  and  full- width  at  half  maximum  (FWHM). 

B.  In  Vivo 

The  MPC,  TTP,  and  FWHM  values  extracted  from  the 
in  vivo  data,  are  listed  in  Table  1  for  all  methods.  The 
methods  were  consistent  (Fig.  4)  in  estimating  the  TTP, 
but  their  results  differed  in  estimating  the  shape  of  HRF 
as  in  simulations.  For  2s  mISI,  the  MPC  of  TW  Averaging 
was  greater  than  the  MPC  of  the  other  two  methods,  which 
is  consistent  with  the  overshooting  of  the  peak  observed  in 
simulations.  For  5s  and  10s  mISI,  the  MPC  of  TW  Av¬ 
eraging  was  less  than  the  MPC  of  the  other  two  methods, 
which  is  also  consistent  with  the  undershooting  of  the  peak 
observed  in  simulations. 

For  the  5s  mISI  data  sets  of  various  randomness,  ex¬ 
tracted  TTP  was  in  agreement  (Fig.  5)  among  the 
three  methods  (mean±SD  =  5.93±0.14,  5.93±0.16,  and 
5.94±0.18).  But  they  differed  in  estimating  the  shape  of 


the  HRF.  Undershooting  or  overshooting  behaviour  of  the 
peak  from  the  TW  Averaging  method  was  similar  to  the 
simulation  results.  Note  that  the  Cepstral  Analysis  did  not 
converge  to  zero  for  type-C  randomness.  The  CG  for  De- 
convolution  method  maybe  more  robust  than  other  meth¬ 
ods  in  that  the  standard  deviations  of  MPC  (0.09)  and 
FWHM  (0.38)  from  the  three  randomness  data  sets  were 
the  smallest  among  the  three  methods. 

IV.  Conclusions 

Our  main  finding  is  that  all  methods  accurately  identified 
the  TTP  of  the  HRF,  but  only  the  CG  for  Deconvolution 
method  was  reliable  in  estimating  the  shape  under  vary¬ 
ing  experimental  conditions.  The  TW  Averaging  was  the 
least  robust  in  estimating  the  shape  of  the  HRF  in  both  the 
simulation  and  in  vivo  data.  Cepstral  Analysis  worked  well 
for  the  zero-noise,  simulated  rER.-fMRI  signals,  however  it 
was  unstable  in  the  presence  of  noise.  The  CG  for  Decon¬ 
volution  method  performed  the  best  both  in  simulations 
and  in  vivo  applications.  Since  it  is  an  iterative  method,  it 
handled  noise  better  than  the  Cepstral  Analysis.  We  con¬ 
cluded  that  the  CG  for  Deconvolution  method  is  a  reliable 
and  robust  method,  and  it  is  a  strong  alternative  to  the 
commonly  used  TW  Averaging  method. 
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TABLE  I 

MPC,  TTP,  and  FWHM  of  Extracted  HRF  from  In  Vivo  Data  for  different  mean-ISI  values 


mean-ISI  |  Cepstral  Analysis  TW  Averaging  CG  for  Deconvolution 


(sec) 

MPC  (%  ) 

TTP  (sec) 

FWHM  (sec) 

MPC  (%  ) 

TTP  (sec) 

FWHM  (sec) 

MPC  (%  ) 

TTP  (sec) 

FWHM  (sec) 

10 

1.23 

6.49 

7.10 

0.88 

6.61 

6.15 

1.22 

7.08 

6.90 

5 

0.83 

6.19 

7.10 

0.77 

6.01 

6.57 

0.79 

6.08 

6.61 

2 

0.44 

6.19 

6.12 

0.49 

6.04 

7.66 

0.43 

5.82 

7.42 

5  type- A 

0.78 

6.00 

7.08 

0.72 

5.90 

6.66 

0.75 

5.96 

6.76 

5  type-B 

0.62 

6.03 

6.15 

0.49 

5.79 

5.91 

0.68 

5.75 

6.06 

5  type-C 

0.51 

5.77 

6.93 

0.69 

6.10 

6.64 

0.58 

6.11 

6.66 

meambSD 

0.64T0.14 

5.93T0.14 

6.72T0.50 

0.63T0.13 

5.93T0.16 

6.40T0.43 

0.67T0.09 

5.94T0.18 

6.49T0.38 

MPC  =  maximum  percentage  change;  TTP  =  time-to-peak;  FWHM  =  full  width  at  half  maximum. 
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Fig.  5.  ( In  vivo  data)  The  effect  of  changing  randomness  on  HRF 

extraction  algorithms.  The  in  vivo  time  courses  for  5s  mISI  are 
shown  on  the  left  column.  Results  of  the  HRF  extraction  algo¬ 
rithms  are  plotted  on  the  right  column  with  the  same  notations 
as  those  in  Fig. 2. 


Appendix 

A.  CG  for  Deconvolution  Algorithm 

Inputs  y:  nxl  signal  vector,  s:  k  x  1  impulse  train  vector. 

Output  h:  m  x  1  least  squares  solution  to  y  =  s  *  h. 
h  =  ones(m,  1)  — >  initialize  h 

h  =  [ h ;  zeros(n  —  m ,  1)],  s  =  [s;  zeros(n  —  k ,  1)] 

Vf  =  fft(v ,  n),  hf  =  fft(h ,  n)  ,  sf  =  fft(s,  n) 

Sf=scfOsf  +  v 

=  sCfQ  yf  -  SfQhf 
r  =  ifft{r_j),  r  =  r(l  :  m),  p  =  — r 


while  ||  r  ||  /  ||  y  ||>  e 

V  =  [p;  zeros(n  -  m,  1)],  pf  =  fft{p,n ) 

— /  =  Sf  OPf  w  =  ifft(wf),  w  =  w(l:m) 

a  =  rTr/pTw 
h  =  h  —  ap 

Told.  =  r,  r  =  rQid  +  aw 
P  =  rTr/r^ldrold 
p  =  —r  +  f3p 
end(while) 

B.  Implementation  Details 

In  CG  for  Deconvolution  algorithm,  we  used  the  follow¬ 
ing  MATLAB  notation,  yf  =  fft(y,n)  is  n-point  FFT  of 
vector  y.  h  =  [h;  zeros(n  —  m,  1)]  means  to  pad  zeros  at 
the  bottom  of  vector  h.  Q  denotes  component-wise  mul¬ 
tiplication.  hc  denotes  component- wise  complex  conjugate 
of  vector  h.  is  represents  the  regularization  parameter. 
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